R Notebook

library(readxl)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables")
path <- read_excel("metabolic_pathways_mastitis_2023.xlsx", sheet = "taxonomy_path")
path_ig <- read_excel("metabolic_pathways_mastitis_2023.xlsx", sheet = "taxonomy_ig")

pathways <- merge(path, path_ig, all.x=T, by=c("Pathway", "Description"))

write.csv(pathways, row.names = F, "/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables/pathways_tax.csv")

##Importing tables

Warning: The working directory was changed to /Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

Aglomerating pathways at the Subclass level

physeq_subclass <- tax_glom(physeq_pathways, taxrank = "Subclass")

Converting phyloseq object to dataframe

write.csv(mean_pathways, "/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/tables/mean_pathways.csv", row.names = F)

converting dataframe to matrix

library(tibble)
path_subclass_mx <- physeq_subclass_df %>% select(sample_ID, Subclass, Abundance) %>% pivot_wider(names_from = Subclass, values_from = Abundance) %>% select(-UNMAPPED,-UNINTEGRATED) %>% column_to_rownames("sample_ID") %>% as.matrix()

Annotations for Heatmap

library(pheatmap)
library(viridis) #color pallet, it's optional
library(RColorBrewer)

#HEATMAP RA USING ln color pallet RBrewer
heatmap <- pheatmap(
  mat               = t(log10(path_subclass_mx+1)),
  border_color      = NA,
  show_colnames     = F,
  show_rownames     = T,
  angle_col = 90,
#  drop_levels       = TRUE,
 # fontsize_col = 4,
 fontsize_row = 5,
#  fontsize          = 14,
 # color             = brewer.pal(9,"RdYlBu"),
#color = inferno(100),  
#number_color = NA,
 annotation_col = annotation_samples,
  annotation_colors = anno_color,
  annotation_names_col = F,
  annotation_names_row = F,
  cluster_cols = T,
  cluster_rows = T,
 clustering_method = "ward.D2",
  gaps_row = FALSE,
)
heatmap

library(ggplot2)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures")
ggsave(plot=heatmap, "heatmap_pathways_mastitis.png", width = 10, height = 5)

Differential abundance analyses

Aggregate phyloseq objects by time points

physeq_day1_week1 <- subset_samples(physeq_subclass, Time_tx%in%c("Day -1","Week 1"))
physeq_day1_week9 <- subset_samples(physeq_subclass, Time_tx%in%c("Day -1","Week 9"))
physeq_day1_week5 <- subset_samples(physeq_subclass, Time_tx%in%c("Day -1","Week 5"))
physeq_week9_5 <- subset_samples(physeq_subclass, Time_tx%in%c("Week 9","Week 1"))
physeq_week9_1 <- subset_samples(physeq_subclass, Time_tx%in%c("Week 9","Week 5"))

Aggregate phyloseq objects by time points

physeq_day1 <- subset_samples(physeq_subclass, Time_tx%in%"Day -1")
physeq_week1 <- subset_samples(physeq_subclass, Time_tx%in%c("Week 1"))
physeq_week5 <- subset_samples(physeq_subclass, Time_tx%in%c("Week 5"))
physeq_week9 <- subset_samples(physeq_subclass, Time_tx%in%c("Week 9"))

Day 0 (other time points were done in the same way)

physeq <- physeq_week5

MaAsLin2

Aglomerating pathways at the Subclass level

physeq_path <- tax_glom(physeq_pathways, taxrank = "Description")

library(microbial)
rf_time <- biomarker(physeq=physeq_path, group = "Time_tx", normalize = F)

Call:
 randomForest(formula = group ~ ., data = data, importance = TRUE,      proximity = TRUE, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 12

        OOB estimate of  error rate: 38.36%
Confusion matrix:
       Day -1 Week 1 Week 5 Week 9 class.error
Day -1     30      7      1      2   0.2500000
Week 1      2     23     12      3   0.4250000
Week 5      1     11     22      6   0.4500000
Week 9      7      1      8     23   0.4102564

rf_time

rf_time %>% head(n=25) %>% ggbarplot(y="Value",x="Description"
                                     #, fill = "Class", palette = "simpsons"
                                     ) + theme(axis.text.x = element_text(angle = 90, hjust = 1))

Converting phyloseq object to dataframe

physeq_df <- psmelt(physeq_path)

top25 <- rf_time %>% head(n=25)
mean <- physeq_df %>% group_by(Description) %>% summarise(Average = mean(Abundance))

order <- physeq_df %>% filter(Description %in% top25$Description) %>% merge(mean, by="Description") %>% mutate(FC=Abundance/Average) %>% group_by(Time_tx, Description) %>% summarise(top = mean(FC)) %>% filter(Time_tx=="Day -1") %>% arrange(top)

top_25_plot <- physeq_df %>% filter(Description %in% top25$Description) %>% merge(mean, by="Description") %>% mutate(FC=Abundance/Average) %>% 
  ggbarplot(x="Description",y="FC", fill="Time_tx", palette = c("#FFD700", "#008000", "#0073C2ff", "#800080"),
            add = c("mean_ci"), position = position_dodge(width = 0.7), orientation="horiz", xlab = "", ylab = "Fold change (mean)",
            #facet.by = "Class",
            order = order$Description)
top_25_plot

library(ggplot2)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures")
ggsave(plot=top_25_plot, "top25_rf_pathways_mastitis.png", width = 10, height = 8)

Treatment plot

physeq_day1 <- subset_samples(physeq_path, Time_tx%in%c("Day -1"))
physeq_week1 <- subset_samples(physeq_path, Time_tx%in%c("Week 1"))
physeq_week5 <- subset_samples(physeq_path, Time_tx%in%c("Week 5"))
physeq_week9 <- subset_samples(physeq_path, Time_tx%in%c("Week 9"))

library(microbial)
rf_treatment <- biomarker(physeq=physeq_week9, group = "Treatment", normalize = F)

Call:
 randomForest(formula = group ~ ., data = data, importance = TRUE,      proximity = TRUE, ntree = ntree) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 3

        OOB estimate of  error rate: 20.51%
Confusion matrix:
           Antibiotic Control class.error
Antibiotic         15       4   0.2105263
Control             4      16   0.2000000

rf_treatment

rf_treatment %>% head(n=25) %>% ggbarplot(y="Value",x="Description"
                                     #, fill = "Class", palette = "simpsons"
                                     ) + theme(axis.text.x = element_text(angle = 90, hjust = 1))

top25 <- rf_treatment %>% head(n=25)
mean <- physeq_df %>% filter(Time_tx=="Week 5") %>% group_by(Description) %>% summarise(Average = mean(Abundance))

order <- physeq_df %>% filter(Description %in% top25$Description, Time_tx=="Week 5") %>% merge(mean, by="Description") %>% mutate(FC=Abundance/Average) %>% group_by(Treatment, Description) %>% summarise(top = mean(FC)) %>% filter(Treatment=="Control") %>% arrange(top)

`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.

top_25_w5_plot <- physeq_df %>% 
  filter(Description %in% top25$Description, Time_tx=="Week 5") %>% 
  merge(mean, by="Description") %>% mutate(FC=Abundance/Average) %>% 
  ggbarplot(x="Description", y="FC", fill="Treatment", palette = "npg",
            add = c("mean_ci"), position = position_dodge(width = 0.7), 
            orientation="horiz", xlab = "", ylab = "Fold change (mean)",
            facet.by = "Time_tx",
            order = order$Description
            )
top_25_w5_plot

physeq_df %>% 
  filter(Description %in% top25$Description, Time_tx=="Week 9") %>% 
  group_by(Description) %>% 
  wilcox_test(Abundance ~ Treatment, paired = F, alternative = "less") %>% 
  add_significance()

physeq_df %>% 
  filter(Description %in% top25$Description, Time_tx=="Week 9") %>% 
  group_by(Description) %>% 
  wilcox_test(Abundance ~ Treatment, paired = F, alternative = "greater") %>% 
  add_significance()

top25 <- rf_treatment %>% head(n=25)
mean <- physeq_df %>% filter(Time_tx=="Week 9") %>% group_by(Description) %>% summarise(Average = mean(Abundance))

order <- physeq_df %>% filter(Description %in% top25$Description, Time_tx=="Week 9") %>% merge(mean, by="Description") %>% mutate(FC=Abundance/Average) %>% group_by(Treatment, Description) %>% summarise(top = mean(FC)) %>% filter(Treatment=="Control") %>% arrange(top)

`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.

top_25_w9_plot <- physeq_df %>% 
  filter(Description %in% top25$Description, Time_tx=="Week 9") %>% 
  merge(mean, by="Description") %>% mutate(FC=Abundance/Average) %>% 
  ggbarplot(x="Description", y="FC", fill="Treatment", palette = "npg",
            add = c("mean_ci"), position = position_dodge(width = 0.7), 
            orientation="horiz", xlab = "", ylab = "Fold change (mean)",
            facet.by = "Time_tx",
            order = order$Description
            )
top_25_w9_plot

top25 <- rf_treatment %>% head(n=25)
mean <- physeq_df %>% filter(Time_tx=="Day -1") %>% group_by(Description) %>% summarise(Average = mean(Abundance))

order <- physeq_df %>% filter(Description %in% top25$Description, Time_tx=="Day -1") %>% merge(mean, by="Description") %>% mutate(FC=Abundance/Average) %>% group_by(Treatment, Description) %>% summarise(top = mean(FC)) %>% filter(Treatment=="Antibiotic") %>% arrange(top)

top_25_d1_plot <- physeq_df %>% 
  filter(Description %in% top25$Description, Time_tx=="Day -1") %>% 
  merge(mean, by="Description") %>% mutate(FC=Abundance/Average) %>% 
  ggbarplot(x="Description", y="FC", fill="Treatment", palette = "npg",
            add = c("mean_ci"), position = position_dodge(width = 0.7), 
            orientation="horiz", xlab = "", ylab = "Fold change (mean)",
            #facet.by = "Class",
            order = order$Description
            )
top_25_d1_plot

w5_w9_diff_plots <- ggarrange(top_25_w5_plot, top_25_w9_plot, labels = c("A","B"), nrow=1, common.legend = T)
w5_w9_diff_plots

library(ggplot2)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures")

Warning: The working directory was changed to /Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

ggsave(plot=w5_w9_diff_plots, "top25_rf_metabolicpaths_mastitis.png", width = 20, height = 8)

Pathways differential

path_diff_plots <- ggarrange(top_25_w1_plot,top_25_w5_plot, top_25_w9_plot, labels = c("A","B","C"), nrow=1, common.legend = T)
path_diff_plots

library(ggplot2)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures")

Warning: The working directory was changed to /Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Metabolomics_mastitis/figures inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

ggsave(plot=path_diff_plots, "top_rf_metabolicpaths_mastitis.png", width = 24, height = 8)